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performed  under  Project  1370,  "Dynamic  Problems  in  Flight  Vehicles,"  Task 
137004,  "Design  Analysis,"  Contract  F33615-75-C-3001.  Capt.  Gerald  Van 
Keuren,  AFFDL/FBR,  was  the  Air  Force  Project  Engineer. 

V.  Y.C.  Young  was  the  principal  investigator  under  the  supervision  of 
M.  R . Brashears. 

The  user's  manual  for  the  computer  program  developed  in  this  study 
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SECTION  I 
INTRODUCTION 

For  aeroelastic  analysis,  the  aerodynamic  forces  are  determined  from 
the  lifting  surface  theory.  In  the  subsonic  regime,  the  theory  is  well  established. 

Most  of  the  approaches  are  based  on  the  kernel  function  method  or  a discretized 
variant  called  the  doublet  lattice  method.  In  contrast,  the  supersonic  lifting  theory 
is  still  under  active  development  and  is  further  complicated  by  the  fact  that  there 
is  a variety  of  possible  methods. 

This  study  is  directed  toward  the  development  of  an  accurate,  efficient,  in- 
expensive numerical  procedure  for  the  solution  of  the  unsteady  pressures  on  a 
planar  wing  undergoing  simple  harmonic  oscillation  in  supersonic  flow.  An  assess- 
ment on  the  theoretical  and  numerical  aspects  of  several  existing  methods  was 
made.  It  was  concluded  that  all  existing  methods  suffer  some  form  of  drawback. 

| 

Accordingly,  an  approach,  which  fully  utilizes  the  merits  of  several  of  the 
various  methods  while  collectively  removing  their  respective  disadvantages,  is 
proposed.  The  method  is  based  on  a finite  element  approach  to  the  supersonic 
kernel  function  method.  It  differs  from  the  supersonic  doublet  lattice  method  in 
the  following  aspects.  Instead  of  trapezoidal  elements,  characteristic  elements 
are  used.  A linear  element  approximation  on  the  lift  is  used  instead  of  the  con- 
stant doublet  strength.  Finally,  the  collocation  point  is  consistently  chosen  to  be 
at  the  nodes. 

*'he  new  method  offers  a unified  theory  for  both  the  subsonic  and  supersonic 
flow.  It  is  also  extendable  to  the  transonic  regime  as  well  as  to  non-planar  appli- 
cations . 
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SECTION  II 
FORMULATION 


2.1  Small  Perturbation  Formulation 

For  a thin  wing  traveling  at  a uniform  velocity,  V,  the  velocity  field 
relative  to  the  Cartesian  coordinates  fixed  to  the  wing  is  given  as 


Vx  = V + u(x,  y,  z,  t) 
Vy  = 0 + v(x,  y,  z,t) 
vz  = 0 + w(x,  y,  z,t) 


where  u,v  and  w are  the  perturbation  velocities  and  satisfy  the  condition 

u/V,  v/V,  w/v  <<  1 

Assuming  the  flow  field  to  be  irrotational,  a perturbation  velocity 
potential,  tj>  , can  be  defined  such  that 


06  0 6 06 

u=s5-  v*37'  w=8l 


According  to  the  small  perturbation  theory,  the  perturbation  velocity 
potential,  <p  , for  the  subsonic,  transonic  and  supersonic  unsteady  flow  is 
governed  by  the  following  general  differential  equation 


1 - M2  - L +6  +6 

V x y xx  yy  Y 


(2  V * f + 0 


where  M is  the  freestream  Mach  number;  y,  the  specific  heat  ratio  of  the 
fluid;  V,  the  freestream  velocity;  and  c,  the  speed  of  sound.  Subscripts 
for  6 denote  differentiation. 
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I 
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This  equation  can  be  further  linearized  for  subsonic  or  supersonic  flow 
(but  not  for  transonic  flow)  to 

(1  - M2)  4 +4>  + <f>  = -V  (2  V <t>  f + O (2) 

xx  yy  zz  2 ' rxt  Ttt'  ' ’ 

Despite  the  similarity  in  form,  the  differential  equation  assumes  quite 
a different  character  depending  on  the  freestream  Mach  number,  M,  For 
M < 1,  the  differential  equation  is  elliptic  while  for  Ivl  > 1,  the  differential 
equation  is  hyperbolic. 


For  the  proposed  study,  Eq.(2)  can  be  rewritten  in  a more  convenient 
form  for  a supersonic  unsteady  flow  as 

(M2  - 1)  <t>  - <t>  - <t>  = - -4  (2  V <t>  f 

' ' rxx  r yy  r zz  c2  xt  Ttt' 

Because  of  its  hyperbolic  nature,  the  equation  is  to  be  satisfied  only  within 
a characteristic  region,  in  contrast  to  the  entire  domain  as  in  the  elliptic 
case.  This  region  is  the  downstream  Mach  cone  originating  from  the  source 
of  distribution.  Conversely,  a receiving  point  in  the  downstream  location 
would  be  affected  only  by  disturbances  generated  within  the  upstream  Mach 
cone  of  the  receiving  point. 


For  an  external  flow,  any  solution  $ to  Eq.  (3)  would  also  have  to  satisfy 
the  boundary  conditions.  For  the  region  ahead  of  the  downstream  Mach  lines 
from  the  leading  edge  of  the  surface  the  freestream  condition  should  prevail. 
The  disturbances  behind  these  downstream  Mach  lines  must  be  directed  out- 
ward from  their  sources.  The  tangency  condition  at  the  surface  is  given  by 


w(x,  y,  0,t)  = jjj:  z(x,  y,  t) 


(4) 


where  w is  the  perturbation  velocity  in  the  z -direction,  and  z(x,  y,t)  repre- 
sents the  instantaneous  displacement  of  the  mean  surface  normal  to  the  x-y 
plane.  D/Dt  is  the  substantial  differential  operator* 


4?  * 


The  pressure  coefficient  is  related  to  <f>  by 


r 2_  D$ 

P~'yzDt 


Equation  (3)  as  written  is  completely  general  within  the  framework  of 
small  perturbation  theory.  It  admits  quite  an  arbitrary  variation  in  time  and 
space  for  <f>  . To  make  the  problem  more  tractable,  a more  restrictive  and 
yet  practical  assumption  of  simple  harmonic  motion  of  the  airfoil  is  usually 
made.  The  condition  is  expressed  as 


z(x,y,t)  = z(x,y)e 


where  i = >F  and  0)  is  the  frequency  of  the  motion.  z(x,  y)  is  called  the  modal 
deflection  which  can  be  prescribed  by  some  polynomial  representation. 


If  the  disturbance  is  due  entirely  to  simple  harmonic  motion,  it  follows 
that  the  effect  as  expressed  by  (ft  would  also  be  harmonic,  i.e., 


(f>(x,y,z,t)  = ?"(x,y,z)  e1 


Substituting  Eq.  (7)  into  Eq.  (3),  one  gets 


<ft  - 0 - <ft 

r xx  yy  rzz 


~J  (2  i(D  V <f>x  - (/<f>) 
c 


where 


(3  = "Vm  - 1 


and  <f>  (x,  y,  z)  is  the  complex  velocity  potential. 


Consistent  with  the  small  perturbation  theory,  the  substantial  differential 
operator  ^ is  given  as: 


■ . *•»,  ' j*.. 


4 


where  second  order  terms  have  been  neglected.  Thus,  the  tangency  condition 
(Eq.  (4))  becomes 


- 


w 


(x.y.o.t)  = v + 


(9) 


and  the  pressure  coefficient  (Eq.  (5))  becomes 


_2_ 

V 


(t + t)  f(x' y ■ 


z) 


(10) 


Equations  (8),  (9)  and  (10)  form  the  framework  for  the  lifting  surface 
theory  formulated  in  terms  t'f  the  velocity  potential. 


2.2  Acceleration  Potential 


From  Euler's  equation,  Prandtl  introduced  the  concept  of  acceleration 
potential 


(11) 


Under  the  small  perturbation  assumption,  4*  also  satisfies  the  same 
linear  differential  equation  as  Eq.(3),  Following  a similar  development  as 
for  <)>,  one  obtains 


' 


If* 


m 


P »F  -if  - vJj 

xx  yy  Tzz 


= — ~t  (2  iu  V J - 


C = 
P 


~T  T U,  y,  a) 


The  use  of  acceleration  potential  (also  often  referred  to  as  pressure  po- 
tential) proves  to  be  especially  simple  in  dealing  with  the  wake  condition  in  the 
lifting  surface  theory. 


Since  no  pressure  discontinuity  can  physically  exist  in  the  flow  field  other 
than  across  the  lifting  surface,  the  acceleration  potential  is  identically  zero  every- 
where except  over  the  planform.  The  complicated  problem  of  accounting  for  the 
wake  can  be  avoided  entirely  by  using  the  acceleration  potential  formulation.  In 
the  case  of  velocity  potential,  a discontinuity  in  velocity  potential  exists  across 
the  wake,  and  the  wake  has  to  be  treated  explicitly.  The  influence  from  the  wake 
is  especially  important  in  the  subsonic  lifting  theory  due  to  the  elliptic  nature  of 
the  flow.  For  supersonic  flow,  with  its  hyperbolic  characteristic,  the  wake  aft 
of  a supersonic  trailing  edge  cannot  exert  any  influence  on  the  lifting  surface,  and 
may  be  neglected  even  if  the  velocity  potential  formulation  is  used.  However,  if 
the  trailing  edge  is  subsonic,  the  forward  Mach  cone  extending  from  a receiving 
point  on  the  trailing  edge  would  include  part  of  the  wake  into  its  region  of  influ- 
ence. Thus,  even  for  supersonic  flow,  the  use  of  acceleration  potential  still  has 
some  definite  advantage  over  the  velocity  potential  -.n  th;  treatment  of  the  wake. 


2.3  Integral  Equation  Approach 


Instead  of  solving  the  linear  partial  differential  equation  directly  in  three- 
dimensional  space  with  the  boundary  conditions  imposed,  the  integral  equation 
method  approaches  the  problem  from  a different  point  of  view.  This  approach 
assumes  that  elementary  solutions  of  some  type  can  be  found,  each  to  satisfy 
the  differential  equation  as  well  as  the  boundary  condition  at  infinity.  Because 
of  the  linearity  of  the  differential  equation,  the  general  solution  can  be  built  up 
from  elementary  solutions  utilizing  the  principle  of  superposition.  Some  typical 
elementary  solutions  are  the  singular  solutions  such  as  source,  sink,  doublet 
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ana  vortex.  The  manner  in  which  these  elementary  solutions  are  combined  is 
dictated  by  the  requirement  that  the  yet  untouched  surface  tangency  condition  be 
satisfied. 


The  larger  the  number  of  these  elementary  solutions  used,  the  better  would 
be  the  general  solution.  Instead  of  taking  a discrete  number  of  these,  we  can  pass 
to  the  continuum  limit  by  replacing  the  summation  over  discrete  quantities  by  the 
integral  over  a spatial  distribution  of  such  quantities.  This  accounts  for  the  ap- 
pearance of  the  integral.  In  summary,  this  can  be  viewed  as  a building  up  of  the 
solution  by  some  known  elementary  solution.  The  problem  of  solving  for  the  sol- 
ution becomes  a problem  of  how  to  combine  the  elementary  solutions  to  achieve 
satisfaction  of  the  surface  boundary  condition. 

The  integral  equation  puts  the  solution  in  an  elegant  form.  An  essentially 
thvee- dimensional  problem  is  reduced  by  one  dimension  to  a surface  integral. 
However,  the  integral  normally  contains  singularities  and  defies  analytical  treat- 
ment. It  is  interesting  to  point  out  that  when  one  resorts  to  numerical  integra- 
tion, one  has  to  essentially  convert  the  continuous  distribution  to  some  discrete 
distribution,  thus  reversing  the  process  used  in  setting  up  the  integral.  The 
whole  lifting  surface  theory  in  the  last  two  decades  was  evolved  around  evalua- 
ting such  an  integral  representation.  Such  an  approach  is  popular  because  while 
the  outer  boundary  extends  to  infinity,  the  quantities  of  interest  are  only  those 
evaluated  on  the  wing  surfaces  such  as  the  surface  pressure  distribution,  etc., 
which  can  readily  be  evaluated  from  the  integral  approach. 


The  integral  equation  can  also  be  derived  more  mathematically  from  the 
differential  equation  by  Green's  function  method.  However,  for  such  a method 
to  be  applicable,  the  differential  equation  must  also  necessarily  be  linear. 

Because  of  the  use  of  superposition  procedures,  the  lifting  theory  is  math- 
ematically limited  to  the  linearized  (small  perturbation)  potential  flow.  As  such, 
it  obviously  has  its  limitation  and  may  not  be  capable  of  fully  describing  the  flow 
field  in  general.  Nevertheless,  due  to  its  simplicity  and  computational  efficiency, 
it  can  serve  its  purpose  quite  well  as  a design  tool  rather  than  a detailed  flow 
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analyzer.  Short  of  using  a three-dimensional  finite  difference  analysis,  which 
would  involve  a prohibitive  amount  of  computer  time,  so  far  there  is  no  better 
nor  more  economical  method  than  the  lifting  theory  in  studying  the  three-dimen- 
sional wing. 
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SECTION  III 
EXISTING  METHODS 

The  object  of  the  lifting  surface  theory  is  to  seek  the  lift  distribution  cor- 
responding to  some  prescribed  downwash  distribution.  The  solution  is  built  up 
from  elementary  solutions  which  satisfy  the  linearized  potential  equation.  For 
the  supersonic  flow,  the  most  commonly  used  building  blocks  are  the  source, 
the  potential  doublet,  the  pressure  doublet  or  a combination  of  such.  Depending 
on  the  type  of  singular  distribution  used,  the  resulting  integral  equation  takes 
on  different  form.  The  existing  methods  for  the  supersonic  lifting  theory  can 
be  categorized  into  four  classes  according  to  the  integral  equation  used. 

3.  1 Integrated  Downwash  Method 

By  superposition  of  pulsating  acoustic  source  distribution  over  the  lifting 
surface  in  a supersonic  flow,  Garrick  and  Rubinow  (Ref.  1)  derived  the  following 
equation 


1 (x> y) = + ^ //w<Vo)k<x-v  y-y0)dxo dv0 


where  K(x-  xq,  y-  yQ)  is  the  kernel  function,  the  detail  form  of  whicn  is  given  in 
Garrick  and  Rubinow  (Ref.  1),  and  will  not  be  repeated  here.  C is  the  region  of 
the  surface  bounded  by  the  forward  Mach  lines  from  the  point  (x,  y). 


As  originally  derived,  Eq,  (14)  is  applicable  only  to  purely  supersonic  flow. 
The  mixed  supersonic  case  with  the  subsonic  leading  edge  must  be  excluded.  To 
extend  the  applicability  to  cover  such  a case,  the  Evvard  concept  of  a diaphragm 
of  unknown  downwash  is  added  to  the  given  wing.  The  combined  wing  and  diaphragm 
region  has  completely  supersonic  edges  and  can  be  treated  by  Eq.(14).  Since 
w(x,  y ) is  not  known  beyond  the  wing  surface,  the  zero  pressure  difference  condition 
instead  of  the  downwash  is  imposed  on  the  diaphragm  region. 


One  of  the  earliest  attempts  to  solve  Eq.  (14)  for  arbitrary  planform  is  the 
"box"  method,  introduced  by  Pines,  Dungundji  and  Neuringer  (Ref.  2).  The  plan- 
form  and  diaphragm  are  divided  into  an  aggregation  of  square  b >xes.  In  each 
box,  the  upwash  is  assumed  to  be  constant  to  facilitate  the  integration  proces s. 

The  method  was  later  improved  by  the  use  of  rectangular  or  rhombic  boxes  with 
sides  along  Mach  lines,  the  so-called  characteristic  boxes.  Another  possibility 
is  the  Mach  boxes  where  rectangular  boxes  with  diagonals  along  Mach  lines  are 
used.  Zartarian  and  Hsu  (Ref.  3)  discussed  the  relative  meri  s and  drawbacks 
of  each  type  of  box.  They  rated  the  Mach  box  as  most  favorable.  In  the  box 
method,  since  the  upwash  is  assumed  constant  over  the  box,  it  can  be  taken  out 
of  the  integral.  The  integration  is  performed  just  on  the  kernel.  This  can  be 
carried  out  over  each  box  in  advance  and  the  results,  called  the  influence  coef- 
ficients, can  be  tabulated  for  a urit,  wq,  as  a function  of  the  position.  Regard- 
less of  what  box  shape  is  used,  difficulties  occur  for  boxes  which  are  adjacent  to  a 
subsonic  leading  edge.  The  upwash  is  singular  in  such  boxes  and  cannot  be  ade- 
quately accounted  for  by  the  constant  value  approximation.  Another  shortcoming 
of  the  method  is  the  jagged  leading  edge  given  by  the  box  approximations,  Moore 
and  Andrew  (Ref.  4)  and  Olsen  (Ref.  5).  This  inexactness  at  the  leading  edge  de- 
finitely affects  flows  downstream  of  the  Mach  cone. 


Since  any  error  at  the  leading  edge  would  propagate  downstream,  efforts 
to  correct  the  jagged  leading  edge  were  reported  by  Woodcock  (Ref.  6).  The  first 
possibility  is  to  take  smaller  boxes  to  reduce  the  jaggedness.  This  improves  the 
accuracy  but  increases  the  computer  time.  Another  approach  is  to  use  a better 
approximation  of  the  upwash  in  the  upstream  part  of  boxes  that  are  cut  by  the  edge. 
The  approximation  is  made  to  possess  the  inverse  square  root  behavior  at  the  lead- 
ing edge.  A third  method  is  to  weight  the  influence  coefficient  by  the  fraction  ofthe 
box  located  on  the  wing. 


In  some  respects,  the  box  method  can  be  classified  as  a primitive  finite  ele- 
ment method.  When  the  shapes  of  the  boxes  are  relaxed  and  the  approximation  of 
the  upwash  within  each  box  is  of  a higher  order  than  constant  values,  we  have  the 
classical  finite  element  method.  Appa  and  Smith  (Ref.  7)  used  the  quadratic 
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non- conforming  triangular  elements  to  calculate  unsteady  supersonic  aerody- 
namic coefficients  from  Eq.(14).  Elements  on  the  planform  need  not  be  the 
same  as  those  on  the  diaphragm.  As  a result  the  elements  can  be  so  arranged 
that  no  partial  element  would  occur  at  the  leading  edge.  The  downwash  is  con- 
tinuous over  the  planform  instead  of  assuming  discontinuous  jumps  as  in  the 
box  method.  All  these  contribute  to  a smooth  pressure  distribution  and  more 
accurate  results.  It  can  thus  be  concluded  that  any  attempt  to  a higher  order 
box"  method  should  be  done  within  the  theory  of  the  finite  element  method. 

The  assessment  of  the  integrated  downwash  method  is  summarized  below. 
Advantages: 

• For  the  purely  supersonic  case,  the  method  is  both  fast  and 
simple.  Instead  of  solving  an  integral  equation,  the  solution 
is  obtained  successively  by  repeated  evaluation  of  the  surface 
integral. 

• By  discretizing  the  planform  into  boxes,  the  method  is  appli- 
cable to  arbitrary  planforms. 

Disadvantages: 


• In  the  box  method,  the  constant  downwash  assumption  was  con- 
ceived for  computational  convenience  on  a desk  calculator.  This 
gives  rise  to  jumps  in  downwash  across  adjacent  boxes,  which 

is  not  physically  meaningful.  Improved  later  by  linear  varia- 
tion in  the  linear  finite  element  approach,  downwash  is  now 
piecewise  continuous  but  still  lacks  smoothness. 

• The  integral  relationship  was  originally  developed  for  the  purely 
supersonic  case.  To  cover  the  mix^d  case,  a diaphragm  must  be 
incorporated  into  the  problem.  The  location  of  this  diaphragm 
region  is  not  unique  for  complex  configurations,  and  its  size 
could  be  almost  as  large  as  the  planform  at  low  supersonic  Mach 
numbers. 


• Since  the  downwash  does  not  vanish  in  the  wake  region,  the  wake 
must  be  explicitly  treated  if  it  is  contained  with  the  forward 
Mach  cone  of  a receiving  point. 

• For  the  mixed  case,  the  velocity  potentials  can  no  longer  be 
solved  successively  but  require  a matrix  solution.  This  is 
because  the  downwash  on  the  diaphragm  is  not  known  and  only 
the  condition  of  zero  lift  prevails  on  the  diaphragm. 
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• Boxes  have  an  inherent  geometric  problem  in  fitting  the  planform. 

The  situation  is  improved  by  the  use  of  finite  elements. 

3.2  Kernel  Function  Method 

Watkins  and  Berman  (Ref,  8)  used  a distribution  of  pressure  doublets  to 
solvt  the  acceleration  potential  equation.  The  resulting  integral  equation  is  of 
the  form 


WU'V)  = o>K<*-V 


where  C is  the  planform  cut  off  by  the  forward  Mach  cone  from  the  point  (x,  y), 
and  K is  the  kernel  function,  the  detail  of  which  is  given  in  Ref.  (8). 

Equatio  i (15)  is  also  applicable  to  the  subsonic  lifting  theory  if  the  super- 
sonic kernel  is  replaced  by  its  subsonic  counterpart.  In  fact,  much  of  the  sub- 
sonic lifting  calculations  are  solved  by  the  kernel  function  method.  Such  a me- 
thod can,  with  some  modification,  be  extended  to  the  supersonic  case.  The  me- 
thod requires  that  the  lift,  I,  be  expressed  as  a linear  combination  of  suitably 
preselected  functions  with  undetermined  coefficients.  The  coefficients  are  de- 
termined by  satisfying  th  a tangency  condition  in  some  approximate  fashion,  e.g., 
collocation.  The  choice  of  the  trial  functions  is  usually  based  on  the  experience 
with  incompressible  flow  and  the  lifting  line  theory.  They  must  vary  in  a cer- 
tain manner  at  the  boundary  of  the  planform  to  satisfy  the  edge  conditions. 

Cunningham  (Ref.  9)  based  his  pressure  series  upon  the  potential  solu- 
tion for  triangular  wings  with  subsonic  leading  edges.  No  separation  of  chord- 
wise  and  spanwise  variables  was  assumed.  Curtis  and  Lingard  (Ref.  10)  chose 
the  pressure  series  similar  to  the  one  in  the  subsonic  case  where  separation  of 
chordwise  and  spanwise  variables  was  assumed.  The  series  is  given  as  a trun- 
cated, double  Fourier  series  multiplied  by  suitable  functions  to  take  into  account 
the  behavior  of  the  pressure  distribution  at  the  leading,  trailing  and  side  edges 
of  the  planform.  Harris  (Ref.  11)  interpolated  the  pressure  series  by  orthogonal 
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polynomials  multiplied  by  some  weighting  function  and  chose  the  interpolating 
points  at  the  zeros  of  the  orthogonal  polynomials.  The  weighting  function  was 
to  assure  proper  behavior  at  edges.  The  interpolation  was  done  under  the  as- 
sumption that  chordwise  and  spanwise  variables  can  be  treated  separately.  Thus, 
it  can  be  said  that  the  choice  is  rather  arbitrary,  which  is  characteristic  of  this 
type  of  semi-analytical  method. 


The  most  costly  operation  is  in  the  evaluation  of  the  surface  integral.  The 
Gaussian- type  quadrature,  modified  to  take  into  account  the  singularity  of  minus 
one-half  power  at  the  leading  edge,  is  usually  employed.  After  the  integration  is 
performed,  the  undetermined  coefficients  are  determined  usually  by  the  colloca- 
tion method.  The  collocation  points  are  usually  determined  on  the  basis  of  two- 
dimensional,  steady  flow  theory.  The  resulting  system  of  algebraic  equations  is 
to  be  solved  for  the  coefficients.  For  N collocation  points,  the  coefficient  matrix 
is  NxN  and  full,  thus  the  solution  becomes  very  expensive  as  N gets  very  large. 


All  the  above  pressure  series  are  smooth  and  are  slow  to  converge  to  the 
true  discontinuous  pressure  distribution.  Recently,  Cunningham  (Ref.  12)  de- 
vised a pressure  series  in  which  the  pressu  e series  based  on  Chebyshev  poly- 
nomials is  weighted  by  some  discontinuous  function.  The  weighting  function  is 
derived  from  supersonic  conical  flow  theory  solution.  This  implicit  discontinu- 
ous characteristic  of  the  pressure  series  gives  a higher  rate  of  convergence 
than  its  conventional  counterpart. 


A variant  of  the  kernel  function  method  is  the  doublet  lattice  method.  The 
planform  is  discretized  with  trapezoidal  panels.  For  subsonic  flow,  a line  of 
constant  doublet  strength  is  assumed  acting  through  the  25%  chord  of  each  panel 
to  approximate  the  lift  distribution.  Collocation  point  is  taken  at  75%  chord  of 
each  panel.  The  subsonic  doublet  lattice  method,  due  to  its  generality,  has  been 
implemented  in  NASTRAN  (Ref.  13).  Harder,  MacNeal  and  Rodden  (Ref.  14) 
showed  that  a straight  forward  extension  of  the  method  to  the  supersonic  flow 
is  not  feasible  due  to  the  singularity  of  the  kernel  along  the  forward  Mach  cone 
of  the  collocating  point.  They  proposed  to  use  Woodward's  constant  pressure 
panel  (Ref.  15)  for  supersonic  application  but  no  result  was  given.  Nevertheless, 
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in  a recent  paper,  Brock  and  Griffin  (Ref.  16)  attempted  the  direct  extension  of 
the  doublet  lattice  method  to  the  supersonic  case.  While  the  collocation  point 
is  kept  at  75%  chord,  it  was  found  that  better  results  were  obtained  for  supersonic 
flow  if  the  line  of  doublet  is  placed  at  50%  chord  of  each  panel.  No  conclusive 
explanation  was  given. 


The  kernel  function  method  is  assessed  below: 


Advantages: 


The  lift  distribution  is  treated  directly  as  the  unknown.  A c‘ 
approximation  to  the  lift  is  sufficient  as  well  as  consistent 
with  the  linearized  supersonic  flow  theory. 


Since  there  is  no  lift  off  the  planform,  it  is  not  necessary  to 
treat  the  wake  nor  to  incorporate  a diaphragm  into  the  analysis. 


The  subsonic  kernel  function  method  is  well  developed.  The 
supersonic  case  retains  essentially  the  same  form.  Since  ker- 
nels have  also  been  formulated  for  transonic  and  nonplanar 
cases,  the  kernel  function  method  offers  a possibility  of  a uni- 
fied approach. 


• The  discretization  as  in  the  doublet  lattice  method  helps  to 
generalize  the  kernel  function  method. 


Disadvantages. 


The  loading  functions  have  to  be  assumed  according  to  the 
planform  and  edge  conditions.  For  any  change  in  planform 
or  in  edge  condition,  a new  set  of  loading  functions  is  re- 
quired. The  method  suffers  from  a lack  of  generality. 


The  loading  functions  used  have  traditionally  been  smooth 
and  continuous,  which  can  hardly  approximate  the  actual 
piecewise  continuous  nature  of  the  loading  in  supersonic 
. w*  situation  can  be  improved  by  incorporating  some 

piecewise  continuous  weighting  function  from  the  conical 
flow  theory. 


A matrix  equation  must  be  solved  to  determine  the  unknown 
coefficients. 


• The  location  for  the  placement  of  the  line  of  doublet  appears 
to  be  arbitrary.  In  any  case,  the  lift  distribution  has  dis- 
continuous jumps  in  values  from  panel  to  panel. 
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3.3  Integrated  Potential  Method 


Based  on  the  generalized  Green' s theorem,  W.  P.  Jones  (Ref.  17)  derived 
an  integral  equation  relating  the  downwash  to  a distribution  of  velocity  potential 
doublets.  Since  the  velocity  potential  doublet  vanishes  on  the  diaphragm,  no  dia- 
phragm  region  needs  to  be  considered.  However,  the  wake  region  still  has  to  be 


treated. 


Allen  and  Sadler  (Ref.  18)  used  a characteristic  mesh  and  approximated 
the  unknown  potential  over  each  rhombus  by  some  linear  or  parabolic  interpo-*- 
lation  involving  the  undetermined  potentials  at  the  nodes  of  the  rhombus.  It  is 
remarkable  that  as  early  as  1963,  they  developed  what  is  now  commonly  known 
as  the  shape  function  in  the  finite  element  method.  Allen  and  Sadler  proceeded 
to  solve  the  nodal  potentials  successively  by  collocation  at  the  nodes,  starting 
at  the  foremost  point  and  progressing  backward  along  the  Mach  lines.  Woodcock 
(Ref.  19)  used  a similar  procedure,  incorporating  minor  improvements. 


Two  recent  integrated  potential  methods  are  those  by  Appa  and  Jones  (Ref. 
20)  and  Giesir-  and  Kalman  (Ref.  21).  Essentially  both  approaches  imply  a lin- 
ear finite  element  approximation  for  the  velocity  potential  doublet  distribution. 
Appa  and  Jones  used  non-uniform  characteristic  mesh  with  the  vertices  as  col- 
location points.  Giesing  and  Kalman  used  trapezoidal  elements  and  set  the  col- 
location point  at  95%  chord,  point  on  the  boxes.  Instead  of  solving  the  unknown 
potential  in  succession,  both  approaches  generate  a set  of  equations  and  the  po- 
tential is  solved  simultaneously. 


The  integrated  potential  method  is  assessed  below. 


Advantages 


• No  diaphragm  is  needed. 

• In  the  approach  by  Allen  and  Sadler,  the  solution  is  obtained 
successively  from  leading  edge  to  trailing  edge. 

• The  formulation  is  more  general  and  is  not  limited  to  planar 
case  alone. 
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Di  sadvantages: 

• The  wake  region  still  has  to  be  treated. 

• The  C°  approximation  of  the  velocity  potential  gives  rise  to 
finite  jumps  or  discontinuity  in  the  pressure  field.  To  obtain 

a meaningful  piecewise  continuous  pressure  field,  a C1  approx* 
imation  of  the  velocity  potential  must  be  employed. 

• In  the  approaches  of  Appa  and  Jones,  and  Giesing  and  Kalman, 
a matrix  equation  must  be  solved. 

3.4  General  Theory  of  Unsteady  Compressible  Potential  Aerodynamics 

Chen,  Suciu  and  Morino  {Ref.  22)  presented  a general  formulation  for  steady 
and  oscillatory,  subsonic  and  supersonic  flow  around  complex  configurations. 
Applying  the  Green  function  on  the  linearized  velocity  potential  equation,  they  ob- 
tained a representation  for  the  potential  at  a control  point  in  terms  of  the  poten- 
tial and  its  normal  derivative  on  the  surface  of  the  body  and  wake.  This  is  equi- 
valent to  a distribution  of  a combination  of  sources  and  potential  doublets.  In- 
stead of  imposing  the  usual  downwash  boundary  condition,  they  obtained  an  inte- 
gral equation  by  requiring  the  potential  to  be  continuous  as  the  control  point  ap- 
proaches the  surface.  A finite  element  approach  with  collocation  is  used  to  ob- 
tain the  solution  by  solving  a set  of  equations. 


This  method  is  assessed  below: 


Advantages: 


• The  formulation  is  unified  for  subsonic  and  supersonic  flow. 

• The  theory  is  applicable  to  non-planar  complex  configurations 
The  thickness  effect  can  be  included. 


Disadvantages: 


• The  wake  and  diaphragm  region  must  be  incorporated  into  the 
analysis. 

• Zeroth  order  finite  element  approximation  is  used  for  the  po- 
tential. e 

9 The  pressure  is  obtained  by  further  differencing  the  potentials 
after  they  are  determined. 

• A matrix  equation  must  be  solved  in  the  solution  process. 


' 


3.5  Recommendations 


Based  on  the  assessments  outlined  in  previous  sections,  it  is  concluded 
that  any  prospective  method  should  have  the  following  desirable  features: 


1.  Eliminate  the  diaphragm  and  wake  region. 

2.  Solve  foi  the  unknown  lift  directly. 

3.  Offer  unified  approach  for  subsonic  and  supersonic  flow,  and 
possibly  transonic  flow. 

) 

4.  Extend  to  non-planar  and  interfering  cases. 

5.  Adopt  finite  element  discretization  to  make  it  applicable  to 
arbitrary  planform. 

6.  Avoid  any  matrix  inversion  from  the  solution  process  to  save 
computer  time  and  storage. 

The  kernel  function  method  comes  closest  to  satisfying  all  the  above  re- 
quirements, except  the  last  two.  In  the  next  section,  the  finite  element  techni- 
que is  applied  to  the  kernel  function  method  to  make  it  more  tractable  by  modern 
computers,  and  help  fulfill  the  last  two  requirements. 


SECTION  IV 
PRESENT  METHOD 

In  this  section,  a new  computational  method  is  presented.  By  using  a 
finite  element  approach  to  the  kernel  function  method,  some  of  the  inherent  short- 
comings of  the  latter  method  are  removed. 

4.  1 Theoretical  Considerations 

The  present  method  is  largely  motivated  by  the  success  of  the  finite  ele- 
ment method  in  replacing  the  Classical  Rayleigh-R  itz  method.  The  difficult  and 
sometimes  impossible  task  of  choosing  the  appropriate  trial  functions  is  conven- 
iently replaced  by  a systematic  and  general  piecewise  polynomial  approximation 
over  the  discretized  domain. 

The  method  is  further  guided  by  the  theory  of  characteristics  which  says 
that  flow  properties,  such  as  the  velocities  and  pressure,  while  continuous  may 
not  have  continuous  first  derivatives  across  characteristic  lines.  Mathematically, 
this  means  flow  properties  are  C°  continuous  in  supersonic  flow.  Since  the  theory 
of  characteristics  is  also  based  on  the  same  linearized  differential  equation,  this 
property  has  to  be  observed  even  in  the  integral  approach.  Thus,  the  lift  distri- 
bution, instead  of  smooth  as  in  the  subsonic  case,  is  now  C°  continuous  for  super- 
sonic flow.  Fortunately,  most  of  the  shape  functions  in  finite  element  method  be- 
long to  the  C class.  In  fact,  it  would  be  much  more  difficult  to  devise  a C1  con- 
tinuous  shape  function, 

| 

Based  on  these  considerations,  the  finite  element  shape  functions  can  ad- 
vantageously be  applied  to  the  kernel  function  method. 

I 
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4.2 


Statement  of  the  Problem 


The  integral  equation  relating  the  complex  upwash  Vw(x,  y)0  to  the  load 
distribution  for  the  supersonic  lifting  surface  theory  is  given  in  the 

AGARD  notation  (Ref.  6)  as 

y)  = 47 '<x0-y0>  K<x- v yy0)dxody0  (16) 

c 

where  C is  the  area  of  the  planform  cut  off  by  the  forward  Mach  cone  from  the 
point  (x,  y). 

The  oscillatory  supersonic  kernel  function  K for  a planar  surface  is  given 
in  the  form 
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Let  the  vertical  displacement  of  the  lifting  surface  be  expressed  as 

Z(x,  y,t)  = s£  f.(x,  y)  q.  (t)  (I8) 


where  s is  the  wing  semi- span,  the  f.(x,  y)  are  the  (dimensionless)  modal  func- 
tions and  the  q.(t)  are  the  (dimensionless)  generalized  coordinates. 

The  nondimensionalized  upwash  on  the  lifting  surface  is 


= (ik+^r)Efi(x’y)qi(t)  <*9) 

i 

th 

Thus  for  the  i mode,  the  corresponding  upwash  is 

wi(x,  y,  0,  t)  = |ik  + ^|-)  f.(x,  y)q.  (t)  (20) 

Substituting  this  in  the  integral  equation,  one  obtains  for  the  ith  mode 

shape 


Equation  (21)  can  be  rearranged  as 


K(x  - x , y - y ) dx  dy 
o 7 7o  o 7o 


(22) 


By  defining  ^(XQ»  yQ)  = ^(XQ»  yo)/<l^(t),  which  can  be  interpreted  as  the 

lift  distribution  due  to  a unit  displacement  in  the  i^  generalized  coordinate, 
the  working  form  of  .e  integral  equation  becomes 


(£  + ik)  fi(x-  y)  yo>  K(x 


x , y - y ) dx  dy 
o 7 7o  o 7o 


(23) 


For  aeroelastic  analysis,  the  aerodynamic  forces  are  conventionally  in- 
put in  the  form  of  some  integrated  force  coefficients  rather  than  actual  lift  dis- 
tribution. The  AGARD  definition  of  the  generalized  force  coefficients  is 


E(x»  y)  Mx,  y)  dx  dy 

J 


(24) 


th  Z 

where  F(x,  y)  are  the  i modal  functions  tnd  pV  Aj(x,  y)  *1^(0  is  the  contribution 
to  the  vertical  force  per  unit  area  resulting  from  harmonic  oscillation  in  the  jth 
degree  of  freedom  (i.e.  with  q.(t)  = q.  Q lwt). 


The  complex  Q. . 

ij 


is  conventionally  written  a;» 


Q..  = Q'..  + ik  Q"  (25) 

y y ij 


where 
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and 

Q1.'.  =Jm(Q..)/k 
ij  ij 

Wltll 

k = us/v,  the  reduced  frequency 

Cjiven  a set  of  modal  functions,  F,  the  geometry  of  the  planform,  the  free- 

stream  Mach  number  and  the  reduced  frequency,  the  problem  is  to  solve  Eq.  (23) 

to  obtain  A^.  This  is  subsequently  input  into  Eq.  (24)  to  compute  the  generalized 

force  coefficients  Q'. . and  Q!1.  . 

lJ  ij 

4.3  Finite  Element  Formulation 

The  first  step  in  the  finite  element  approach  is  to  discretize  the  planform 
into  finite  elements.  With  arbitrary  aligned  elements,  since  the  region  of  inte- 
gration is  limited  to  the  forward  Mach  cone,  elements  on  the  Mach  cone  boundary 
would  be  partially  cut  by  the  Mach  cone  (Fig.  1).  The  partially  cut  elements  would 
have  to  be  determined  prior  to  the  integration  process  for  each  collocation  point. 
Clearly  this  would  involve  a considerable  amount  of  bookkeeping,  not  to  mention 
the  labor  in  determining  the  intersections.  On  the  other  hand  if  all  the  elements 
are  bounded  by  characteristic  lines,  it  becomes  a clearcut  decision  as  to  whether 
the  element  is  to  be  included  in  the  integration.  For  this  method,  a characteris- 
tic mesh  would  be  used  exclusively. 


After  the  planform  is  discretized  into  elements,  the  surface  integral  can 
be  replaced  by  a summation  of  surface  integrals  over  each  element  contained 
within  the  forward  Mach  cone.  Equation  (23)  becomes 


where  i8  the  local  lift  distribution  over  the  element  and^  denotes  summa- 
tion over  elements  within  the  forward  Mach  cone.  The  local  lift  distribution  can 
be  replaced  by  the  finite  element  approximation  as 

Xie,(xo-yo)  =fTt(xo.yo).A[e)  (27) 

where  N is  the  row  vector  of  shape  function  and  is  the  column  vector  of 
nodal  lift  values. 

Substituting  Eq.  (27)  into  Eq,  (26),  one  obtains  the  finite  element  form 
of  the  integral  equation 

+ ik)  fi(x*  = iX ff  n 1 <v  yQ]  K(x  ■ v y ■ y0)dA  • *le)  <28> 


One  immediate  consequence  is  that  the  integrand  no  longer  contains  the  un- 
known lift  distribution.  It  depends  only  on  the  relative  position  between  the  influ- 
encing element  and  the  collocation  point,  and  can  be  evaluated  readily  with  pro- 
per attention  to  the  singularities  in  the  kernel  function. 


For  convenience,  define  the  integrated  values  of  the  product  of  shape  functit 
and  kernel  function  as  the  weighted  kernel  coefficients,  i.e.. 


c t<e)  ff  <xo*  yo)  K(x ' v y ' y0)dA 


The  integral  equation  finally  becomes 


(a7  + ik)  V*.yl  = 
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With  the  finite  element  approximation,  the  generalized  force  coefficients 


are 


N 


Q ..  = - 

ij 


T,  ff  E(x,  y)  (x,  y)dxdy  . 

n = \JJ(e)  J 


(31) 


N 


where ^ denotes  summation  over  all  the  elements  on  the  planform. 
n - 1 


Equations  (29),  (30),  and  (31)  represent  the  finite  element  formulation 
of  the  lifting  surface  theory. 


4.4  Choice  of  Elements 


Originally  it  was  thought  that  by  using  higher  order  triangular  and  quadri- 
lateral elements,  fewer  nodes  are  required  to  generate  results  of  the  same  ac- 
curacy. While  this  is  true  in  general,  it  is  not  suitable  for  the  type  of  solution 
process  used  in  this  method. 


Since  the  nodes  are  to  be  used  as  collocation  points,  the  additional  nodes  on 
the  side  would  have  their  forward  Mach  cones  cutting  through  the  elements  (Fig.  2), 
Also,  since  the  nodal  variables  of  the  same  element  are  dependent  on  each  other, 
all  the  unknown  nodal  variables  would  have  to  be  solved  simultaneously  for  each 
element.  For  example,  if  quadratic  elements  are  used,  at  each  stage  there  would 
be  three  unknown  nodal  variables  to  he  determined  simultaneously.  Whereas,  if 
linear  elements  are  used,  there  would  be  just  one  unknown  nodal  variable  to  be 
determined  at  each  stage.  The  net  saving  by  using  higher  order  elements  in  elim- 
inating the  interior  nodes  can  hardly  justify  the  additional  work.  Also,  the  para- 
metric property  of  the  higher  order  elements,  which  can  be  used  to  fit  curve  char- 
acteristics, can  hardly  be  utilized  in  this  study  since  the  Mach  number  remains 
constant  over  a supersonic  planar  surface.  With  due  consideration  to  all  these 
problems,  it  was  concluded  that  the  linear  finite  elements  can  best  serve  our  sol- 
ution technique. 
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4.5  Choice  of  Mesh 


The  nonuniform  element  mesh  was  initially  considered  since  it  would 
provide  an  added  flexibility  in  fitting  the  actual  planform  (Fig.  3).  While  this 
is  technically  feasible,  it  was  not  adopted  due  to  economic  reasons.  For  each 
collocation  point,  every  element  contained  in  the  forward  Mach  cone  must  be 
integrated.  The  surface  integral  is  usually  done  numerically  through  the 
use  of  Gpussian  quadratures.  Since  the  integrand  contains  the  kernel  function 
with  its  singularities,  a high  order  Gaussian  quadrature  has  to  be  employed 
to  get  a meaningful  representation.  As  one  marches  downstream  to  the  trail- 
ing edge,  more  and  more  elements  are  included  in  the  Mach  cone.  The  costly 
numerical  integration  has  to  be  performed  over  and  over  even  for  the  same 
element  since  the  relative  position  of  the  receiving  point  is  changing  each 
time.  For  a coarse  mesh,  this  does  not  present  much  of  a problem.  If  the 
mesh  has  to  be  refined  for  better  accuracy,  which  is  almost  mandatory  for 
unsteady  flow  at  high  frequency,  the  wastefulness  of  such  nonuniform  mesh 
becomes  clear.  In  an  effort  to  reduce  the  computer  time,  it  was  realized 
that  since  the  kernel  function  is  a function  only  of  the  relative  location  between 
the  receiving  point  and  the  sending  point,  much  of  the  duplicated  effort  in  the 
numerical  integration  can  be  avoided  by  using  a uniform  mesh.  Essentially 
a stencil  of  uniform  elements  large  enough  to  cover  the  most  extreme  cases 
is  set  up  and  the  integrated  values  are  stored  in  a table.  For  each  receiving 
point,  the  integration  over  regular  elements  can  be  replaced  simply  by  a table 
look-up,  with  the  necessary  numerical  integration  performed  only  on  the  irreg- 
ular elements.  With  such  an  approach,  much  saving  in  computer  time  can  be 
realized. 


For  the  rectangular  wing,  since  the  planform  edges  are  composed  of 
either  horizontal  or  vertical  lines,  the  mesh  can  be  arranged  such  that  the 
edges  cut  the  characteristic  boxes  regularly  across  the  vertices.  As  a result, 
only  a combination  of  characteristic  boxes  and  regular  triangular  elements 
is  necessary  to  fit  the  planform  exactly  (Fig.  4).  However,  for  the  non- 
rectangular  planform,  since  the  edges  are  inclined  at  arbitrary  angles,  they 
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cut  the  boxes  across  the  sides  at  several  possible  locations  (Fig.  5).  This 
gives  rise  to  the  possibility  of  having  a partially  cut  element  with  five  vertices. 

Some  effort  was  spent  in  developing  subroutines  to  handle  these  irregu- 
lar polygons.  If  the  inverse  square  singularity  does  not  pass  through  these 
elements,  they  can  easily  be  treated.  For  example,  the  five-sided  element 
can  be  treated  as  a combination  of  triangle  and  quadrilateral.  However,  for 
singular  elements,  some  additional  problems  arose.  Since  the  width  of  the 
singular  strip  is  limited  by  the  proximity  of  the  nodes,  some  situation  might 
arise  where  the  strip  has  to  be  taken  extremely  small.  Tests  on  the  sub- 
routines indicated  that  for  such  extreme  cases,  the  nodal  values  of  the  inte- 
grated kernel  can  become  extremely  large  in  magnitude  and  unreliable.  It 
should  be  pointed  out  here  that  in  Appa's  finite  element  approach  on  the  inte- 
grated downwash  method  (Ref.  7),  it  was  possible  to  fit  the  planform  exactly 
mainly  because  the  only  singularity  is  of  the  inverse  square  root  type,  which 
is  integrable,  and  not  of  the  inverse  square  type,  which  necessitates  the  use 
of  Cauchy  principal  value. 


To  keep  the  program  simple,  it  was  then  decided  to  use  just  the  type  of 
elements  which  have  been  proved  successful  with  the  rectangular  wing.  For 
the  non-rectangu'ar  wing,  this  gives  rise  to  the  jagged  leading  edge,  a common 
undesirable  feature  of  all  box  methods.  However,  with  the  finite  element  tech- 
nique, by  introducing  triangular  fill-in  elements,  the  jaggedness  at  the  planform 
edge  is  substantially  reduced.  The  method,  in  effect,  combines  the  best  features 
of  the  Mach  box  method  and  the  characteristic  box  method.  For  example,  it 
can  approximate  a rectangular  wing  as  exactly  as  the  Mach  box  method,  while 
it  can  also  handle  sonic  edges  as  well  as  the  characteristic  box  method. 


4.6  Finite  Element  Integration  Scheme 

lim 

The  supersonic  kernel  function  has  a singularity  of  the  type  e 

e 0 

along  the  line  extending  from  the  collocation  point  in  the  upstream  direction. 


I 


This  inverse  square  singularity  is  not  integrable  in  the  usual  sense,  and  the  im- 
proper integral  must  be  evaluated  using  the  concept  of  Cauchy's  principal  value. 

i 1 

tv  im  \T" 

ihe  supersonic  kernel  function  has  an  additional  singularity  of  the  type  v<r 

along  the  forward  Mach  cone  boundary.  Fortunately,  this  type  of  singularity  is 
integrable  and  can  easily  be  treated  by  the  Gaussian  type  quadrature. 


In  the  conventional  kernel  function  method,  the  area  of  integration  is  divided 
into  sub- regions  such  that  the  spanwise  singularity  is  confined  within  a narrow 
strip  (Fig.  6).  For  regions  other  than  this  strip  the  integral  is  not  singular  and 
the  integration  can  be  carried  out  individually.  Care  must  be  taken  to  apply  the 
proper  integration  limits  which  are  dependent  on  the  planform  geometry  and  the 
location  of  the  control  point.  Treatment  of  the  singular  strip  falls  into  two  cate- 
gories. The  first  method  is  to  carry  out  the  chordwise  integration  in  conjunction 
with  a Lagrange  interpolation  in  the  spanwise  direction.  The  resulting  polynomial 
is  then  integrated  using  the  Cauchy  principal  value  to  extract  the  finite  part  of  the 
improper  integral.  The  second  method  is  to  subtract  and  add  a term  which  mimics 
the  integrand  at  the  singularity.  The  original  integrand  with  the  singularity  sub- 
tracted is  now  well  behaved  and  can  be  evaluated  numerically.  The  correction  term 
which  contains  a simpler  integrand  can  be  evaluated  analytically  using  the  Cauchy 
principal  value. 


In  the  finite  element  approach,  since  the  local  approximation  is  used,  it 
would  be  more  efficient  to  perform  the  integration  on  an  element- to- element  basis. 
Since  a uniform  characteristic  mesh  of  linear  elements  is  used,  the  nodes  are  ar- 
ranged in  a regular  pattern.  The  singular  strip  is  taken  to  be  a fraction  of  the  ele- 
ment half  width.  If  the  singular  strip  passes  through  an  element,  it  does  so  regu- 
larly through  the  middle  nodes  or  through  either  extreme  nodes.  (Fig.  7).  For  such 
singular  elements,  the  integration  is  performed  by  first  decomposing  the  element 
into  sub-regions  and  then  applying  the  appropriate  integration  scheme.  For  non- 
singular elements,  a direct  Gaussian  integration  is  used 
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4.7  Integration  Across  the  Singular  Strip 


t 


As  mentioned  in  the  last  section,  two  methods  are  available  for  the  integra- 
tion across  the  singular  strip.  Due  to  its  simplicity  in  implementation,  the  first 
method  is  adopted. 

The  general  form  of  the  integral  across  the  singular  strip  is 


rVy+€  /.yn) 

j J 

r\3L=yme  ^a(7l) 
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(32) 


K is  the  kernel  function,  which  can  be  written  as 


K = 
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(33) 


(y-nr 


where  the  inverse  square  singularity  has  been  factored  out  and  K contains  only 
integrable  singularities.  Equation  (32)  can  be  written  as 


f , 

f, 

y-f  l-4(n) 


i'l1  (|»n)  K(x,  y; ^ , n , u,M)d£  dq 


(34) 


The  integral  can  be  evaluated  first  b^  .ntegrating  in  the  chordwise  direc- 
tion and  then  in  the  spanwise  direction  q. 

The  chordwise  integration  can  be  replaced  by  some  Gaussian  type  quadrature 
in  the  following  manner 
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Cunningham  (Ref.  9)  presented  a tenth  degree  version  of  the  following 


/ ^7 

J <y-n> 


= 57T33S7  [33’?“<F1  + FU>  + 266.50°IF2  + F10) 


+ 147,375(F3  + F9)  + 822.000(F4  t Fg) 


+ 4,562,250(F5  + F?) 


4 (-12,807,144)  F, 


The  abscissas  are  numbered  starting  from  (y-e  ) to  (y-K  ).  F.,F, , . . . , 

1 u 

Fr  represent  the  function  F (rj ) evaluated  at  the  corresponding  abscissas. 


Tests  showed  that  the  sixth  degree  quadrature  performs  nearly  as  well  as 
the  tenth  degree  version.  For  efficiency,  the  sixth  degree  quadrature,  Eq.  (36), 
is  used  in  present  study. 


4.8  Starting  Solution 


For  a hyperbolic  differential  equation,  the  initial  condition  is  required  to 
start  the  solution  process,  as  in  the  method  of  characteristics.  For  the  integral 
approach,  a starting  solution  is  seldom  used.  In  order  to  take  advantage  of  the 
hyperbolic  nature  of  the  flow  and  also  to  avoid  any  matrix  inversion,  a successive 
solution  process  with  a starting  solution  is  proposed. 


Since  the  dependent  variable  of  the  kernel  function  method  is  the  lift,  the 
starting  olution  requires  that  the  lift  at  the  leading  edges  be  specified.  Adja- 
cent to  the  supersonic  leading  edge  the  flow  is  locally  two-dimension?.!  for  an 


V — T 


oscillating  airfoil  in  supersonic  flow,  Miles  (Ref.  24) 
tude  as 


gave  the  pressure  ampli- 


'W-a0gMv(0+)+a0jf  eU-S)[v5(|)tikv(|)]dg  (38) 

*h"'  30  l 4/3  " 'h,e  tW°‘dirr'"'a,0"al  (Ackeret)  lift  curve  slope,  v i.  the  dimen- 
s ‘unless  downwash  (positive  down),  and 


g(x)  = exp(-ikM2x/p2)J0(kMx/(32)  {39) 

where  JQ  is  the  Bessel  function. 

In  AOARD  convention,  lift  is  normalised  b„  pV2  instead  of  the  usual  pv2/Z 

lo  account  for  the  difference  a - ? /a  ■ , _ . ' 

“erence.  aQ  - 2/0  is  used.  Setting  x = 0,  the  lift  at  a super- 
sonic leading  edge  is  obtained  as  F 
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Replacing  the  downwash  by  our  upwash.w,  Eq. 


(40)  becomes 
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For  a supersonic  swept-back  wing  of  infinite 


a sweepback  factor  (Ref.  25) 


span,  the  lift  is  modified  by 


2 - 1 

~s  w . i 


with  n = A 

0 

where  A is  the  sweepback  angle. 


w "V  * — 77 
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Some 


1 ' 1 
lim  — 

For  a subsonic:  leading  edge,  the  lift  goes  to  infinity  as  \K 

arbitrary  large  number  can  be  specified.  In  this  case,  it  is  very  important  to 
use  the  appropriate  Gaussian  type  quadrature  with  the  same  singular  behavior 
when  performing  the  integration  over  the  leading  edge  element. 


Edge  Conditions  and  Gaussian  Quadrature 


In  arriving  at  the  integral  equation,  the  boundary  condition  at  infinity  is 
satisfied  by  the  judicial  choice  of  the  elementary  solution  such  as  the  source  or 
doublet.  The  nature  of  the  equation,  be  it  elliptic  or  hyperbolic,  is  taken  into  ac- 
count by  observing  the  proper  region  of  integration.  The  tangency  condition  at 
the  surface  is  explicitly  expressed  in  terms  of  the  integral  equation.  However, 
no  provision  was  made  for  the  edge  condition.  This  latter  condition  must  be  in- 
corporated as  a constraint  to  the  integral  equation. 


Edges  are  designated  as  subsonic  or  supersonic  edges  according  to  whether 
the  normal  component  of  the  freestream  velocity  to  the  edge  is  subsonic  or  super- 
sonic. The  edge  conditions  for  Ap  are 

(T)  lim  \/e~  at  side  edge,  subsonic  or  sonic 
”0  trailing  edge 

©<i  m 1 A/e  at  subsonic  leading  edge 

®AP  finite  at  supersonic  leading  edge  or  super- 

sonic trailing  edge. 


These  conditions  can  be  incorporated  into  the  solution  in  the  form  of  weight- 
ing functions  to  the  edge  elements. 
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The 


fim  1/  Ve~  . , , 

c —AO  ^Pe  of  singularity  can  be  handled  readily  by  using  the 
Gauss-Chebyshev  quadrature  as  follows 


where 


= w.f(x.) 

i=l 


X. 

1 


-cos 


t2i  - ill 


(43) 


w . 

1 


1L 

n 


When  applied  to  a function  F(x>  which  contain,  implicit  singularity  of  type 


1 


= 2 wiV^^F(x.) 

i = l 

= £ wiF<*i> 

i = l 
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V*l- 


M 


i 


where 


(2i - 1 ) n 


W.  = w.  VI  -x^ 


[2i  - l)y 


One  of  the  nicest  features  of  the  Gauss -Chebyshev  quadrature  is  that  the 
weights  and  abscissas  can  be  generated  systematically  to  any  degree,  instead  of 
requiring  a table  input  as  with  the  rest  of  the  Gaussian  quadrature.  In  addition, 
it  has  been  determined  in  this  study  that  the  modified  Gauss-Chebyshev  quadrature 
can  also  adequately  approximate  integrals  of  the  type 


J f(x)  - x 


f(x)  dx 


For  f(x)  - 1,  x or  x , less  than  1%  relative  error  can  be  achieved  with  n > 12. 

For  ease  of  implementation,  it  was  decided  to  use  this  "all  purpose"  qua- 
drature, Eq.  (44),  throughout  despite  the  fact  that  relatively  large  numbers  of 
quadrature  points  are  necessary. 

In  the  solution  process,  the  nodal  lift  values  are  initially  set  to  zero.  Nodes 
on  the  side  edges  are  excluded  from  further  treatment  and  retain  their  zero  lift 
value^oughout  the  analysis.  At  the  leading  edge,  the  sweepback  factor, 

1/Vl  - n is  first  determined.  For  0 < n < 1,  the  leading  edge  is  supersonic 
and  the  lift  at  the  leading  edge  is  given  by  Eq.  (42),  where  the  unswept  lift  is 
amplified  by  the  sweepback  factor.  For  n > 1,  the  leading  edge  is  subsonic  and 
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the  sweephack  factor  becomes  imaginary.  The  lift  becomes  singular  at  the 

subsonic  leading  edge  as  lim  C/%  e , where  C is  an  arbitrary  constant  of  0(1). 

e — - 0 

Comparing  this  form  with  the  expression  of  Eq.  (42),  one  can  let  the  sweep- 

back  factor  be  replaced  by  the  factor  lim  1/Ve"  for  n > 1.  It  follows  that 

e — 0 

C is  given  by  the  supersonic  unswept  value  as  expressed  by  Eq.(41).  Thus, 
for  a subsonic  leading  edge,  the  starting  value  is  given  by  Eq.  (41)  while  the 
factor  l/Ve-  is  absorbed  into  the  quadrature  as  a weighting  function. 

For  a supersonic  trailing  edge,  the  computation  is  performed  up  to  and 
including  the  trailing  edge  nodes.  For  subsonic  or  sonic  trailing  edge,  the 
trailing  edge  nodes  are  excluded  from  the  analysis  and  thus  retain  the  zero 
lift  value.  This  means  that  the  lift  distribution  drops  from  a finite  value  to 
zero  within  the  trailing  edge  element. 

4.10  Kernel  Function  Evaluation 

The  kernel  function  expression  in  Eq.  (17)  involves  an  integral  which 
cannot  be  evaluated  analytically.  With  some  numerical  approximation  on 
this  integral,  closed  form  evaluation  of  the  kernel  function  is  possible.  H.  J. 
Cunningham  presented  a closed  form  evaluation  of  Eq.  (17)  in  the  appendix 
of  Ref.  (26). 

Harder  and  Rodden  (Ref.  27)  obtained  a more  general  form  of  the  kernel 
function,  which  is  reducible  to  Eq.  (17)  for  the  planar  case.  A closed  form 
evaluation  of  the  version  by  Harder  and  Rodden  was  given  by  A.M.  Cunningham 
in  the  appendix  of  Ref.  (12). 

In  the  course  of  this  study,  it  was  found  that  some  discrepancies  exist 
in  the  values  of  the  kernel  function  obtained  by  either  procedure.  Some  repre- 
sentative values  are  presented  in  Table  1,  which  shows  that  A.M.  Cunningham's 
version  is  less  sensitive  to  changes  in  reduced  frequency.  Both  versions  were 
subsequently  applied  to  a rectangular  wing  at  M = 1.2.  Results  based  on  A.M. 
Cunningham's  version  gave  better  agreement  to  available  data  (Ref.  6), 


■ , 


especially  for  nonzero  reduced  frequencies.  With  the  added  advantage  of  e 
tensibility  to  nonplanar  case,  A.M.  Cunningham's  version  is  adopted  in  this 
study. 


4.11  Solution  Procedure 


The  solution  procedure  of  this  method  is  summarized  as  follows 


1.  Based  on  the  Mach  number  of  the  flow,  set  up  and  input  the  element 
information  for  the  planform. 

2.  Input  a reduced  frequency  and  the  shape  modes  fi  of  the  deflection 
to  compute  the  nodal  upwash  for  each  mode. 

3.  Clear  nodal  lift  values  and  impose  starting  solution  at  leading  edge 
and  set  other  side  conditions. 

4.  Select  the  collocation  node  (NODE)  starting  with  the  foremost  node 
with  unknown  lift.  Set  SUM  i = 0. 

5.  Select  the  influencing  element. 

6.  Compute  the  weighted  kernel  coefficients. 

7.  SUM.  = SUM.  + C t(e)  . Xje) 

8.  Repeat  steps  5 through  7 until  all  influencing  elements  have  been 
accounted  for. 

9.  Compute  the  unknown  lift  for  each  mode,  i,  according  to 

A.(NODE)  = (4ttw.  - SUM .) /c(NODE) 

10.  Repeat  steps  4 through  9 until  all  nodes  with  unknown  lift  are 
accounted  for, 

1 1.  Compute  the  generalized  force  coefficients. 

12.  Output. 

13.  End. 
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SECTION  V 


RESULTS  AND  DISCUSSION 


A computer  program  (Ref.  28)  based  on  the  present  method  was  developed 
and  applied  to  the  AGARD  planforms  which  include  the  rectangular  A = 2.0,  the 
arrowhead  A = 4.0,  and  the  tapered  swept  back  A = 1.45.  The  choice  of  Mach 
number  and  planform  is  such  that  all  possible  combination  leading/trailing 
edge  conditions  are  covered  (see  Table  2).  The  program  was  run  on  a UNIVAC 
1108  and  the  statistics  given  in  Tables  3 through  5 are  based  on  this  machine. 
The  15  point  modified  Gauss -Chebyshev  quadratured  was  used  throughout  ex- 
cept for  the  Tapered  Sweptback  A = 1.45,  M = 1.04  case,  where  only  12  points 
were  used.  In  general,  a factor  of  2.5  to  3.0  reducation  in  time  can  be  expected 
with  a CDC  6600,  with  a slight  increase  in  accuracy  due  to  the  larger  word  size 
in  reducing  the  truncation  error  in  accumulation. 


Results  are  given  in  the  form  of  generalized  force  coefficients.  The 
most  extensive  source  of  generalized  force  coefficients  is  the  AGARD  report 
by  Woodcock  (Ref.  6).  With  the  exception  of  a few  cases,  all  the  supersonic 
data  in  the  report  were  generated  from  the  characteristic  box  method  (Refs, 
^9  and  30).  Another  source  of  data,  based  on  the  Mach  box  method,  was  pro- 
vided in  the  report  by  Olsen  (Ref.  5).  No  results  based  on  the  kernel  function 
method  are  available  for  comparison.  It  appears  that  this  work  is  the  only 
one  in  tabulating  the  generalized  force  coefficients  by  the  supersonic  kernel 
function  method.  Thus,  the  two  reports  mentioned  earlier  will  serve  as  the 
standard  of  comparison  for  this  study. 


The  generalized  force  coefficients  QL  and  Q'L  are  tabulated  individually 
for  each  case  preceded  by  a figure  depicting  the  actual  mesh  used  for  each 


Mach  number.  One  case  is  defined  as  one  planform,  at  one  Mach  number  end 
reduced  frequency.  Instead  of  separating  the  modes  into  symmetric  sets  and 
antisymmetric  sets  as  in  Rei.  6,  only  one  combined  set  is  used.  It  is  pointed 
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out  that  in  Ref.  5 the  program  used  could  not  handle  antisymmetric  modes  y 
and  x.  Consequently,  the  antisymmetric  modes  were  replaced  by  |y|  and  x |yj 
in  Ref.  5.  One  should  take  this  into  consideration  when  comparing  the  results. 

Also,  the  flap  rotation  mode  in  Ref.  6 is  not  treated  in  this  study. 

Each  table  is  arranged  such  that  CR . is  tabulated  before  Q'.'.  . The  plan- 
form,  Mach  number,  reduced  frequency  and  mode  shapes  are  als^  specified. 

In  addition,  the  mesh  fineness  is  also  specified  through  the  length  of  the  side 
of  the  characteristic  element,  DELTA. 

Because  of  the  large  amount  of  data  involved,  it  is  not  possible  to 
discuss  each  case  individually.  In  general,  the  present  results  agree  reason- 
ably well  in  trends  with  reported  values.  In  some  cases,  the  values  fall  be- 
tween the  reported  values,  while  in  other  cases,  the  values  are  slightly  off. 
Considering  that  discrepancies  exist  even  between  the  two  versions  of  the 

characteristic  box  method,  the  present  results  are  very  encouraging. 

i 

Some  of  the  discrepancies  can  be  explained  readily  from  a mathematical 
point  of  view.  The  generalized  force  coefficients  can  also  be  interpreted  as 
the  moments  of  the  lift  distribution.  In  the  box  method,  the  lift  distribution  is 
approximated  by  some  step  function  like  distribution,  since  the  lift  is  constant 
within  each  box  but  assumes  different  values  from  one  to  the  other.  For  this 
type  of  distribution,  the  approximation  is  accurate  at  most  to  the  zeroth  moment. 
Higner  moments  will  become  less  and  less  accurate.  In  comparison,  the  linear 
approximation  in  the  present  method  is  accurate  up  to  the  first  moment. 

In  conclusion,  since  all  the  reported  values  are  based  on  some  type  of 
box  method,  they  tend  to  reinforce  each  other.  Without  further  unbiased  com- 
parison from  some  other  method,  it  is  almost  impossible  to  conclude  which 
method  gives  the  better  result. 

i 
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Regular  Characteristic  Element 


Triangular  Fill-In  Elements 


Fig.  4 - Typical  Elements  Used  to  Fit  the  Planform 
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Fig.  5 - Regular  Mesh  with  Irregular  Elements  at  Planform  Edges 
(Note  the  occurrence  of  elements  with  five  vertices). 
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Table  1 

COMPARISON  OF  THE  TWO  VERSIONS  OF 
KERNEL  SUBROUTINE 


M = 1.2 

x = O.J 
o 

Vq  s 0,1 


Reduced 
F requency, 
k 

Modified  Kernel, 

K = y2  K/2 
o 

Herbert  Cunningham 
(Ref.  26) 

Atlee  Cunningham 
(Ref.  12) 

0.0 

(-1.3363,  0.) 

(-1.3363,  0.) 

0.1 

(-1  .3359,  .026725) 

(-1.3361,  .019948) 

1.0 

(-1.2969.  .26184) 

(-1.31  16,  .19596) 

45 


rv " ■' 


re-re  ^ . 


G 

VH  O . 

0 iK 
di  flj  * 

O *H 

(J  -H 

^ Ifl  **-' 

O 


(M 

5k 

CM 

X 

«k 

•L 

CM 

(M 

3k 

3k 

•> 

* 

CM 

(M 

X 

X 

* ^ 

» 

X X 

X X 

•>  •» 

— * >S 

**  >* 

H 

•» 

vO 

• 

•» 

CO 

ITl 

• 

• 

«k 

% 

«>  n 

u * u 

« 9 -H 

CL,  CX  c 

3 3 0 

WWW 


u u u 

1 I 1 „ 

O O O .S 

CA  CO  (A  C 


*<  3<  o 

ii  «j  II  n 

a a a 43 


3 3 3 3 

W W W W 


a>  <u 

a a 

3 3 

w w 


G G o 

o o .54 

«0  «5  G 

Lh  Ln  O 

V V CA 

a a xi 

3 3 3 

WWW 


u o u S 

H .H  .H  O 


O O h o 

CA  CO  ft)  i 

•o  £>  a 43 

3 3 3 3 

w w w w 


u 

u 1 

•H  O 

G W 

O h 

CA  ft) 

43  a 

3 3 

w w 


CM 

cm  in  vo 

N O H N Ifl  O O 

••••••• 

^ N H H H N H 


V 

1 " 

2 < 

h 

< 


a 

V 

I !} 
1 

14  II 

V . 

a < 

flj 

H 


RECTANGULAR  A*2  * • N* 1 , 2 • K*.J,  DELTA*. Q7S,  NUDES* 


2.7999-08  1.9S7S-07  6.3058-10  -|.0l22-07  -3.2918-08  -|.923l-0t  2.9529-Ql 


RECTAN6ULAR  A-2.,  N-1.2,  Kb) 


i 


MnnM 

MM0MHW 


fr.4v.,**.*  5K 


000 


RECTANGULAR  A»2..  M*2..  <*.3,  DEL T A • • 0833  . MODES*  1.  X.  X»X,  Y *Y,  X»X«Y«Y 


RECTANGULAR  «>2>|  M"2.,  *■.*,  DELTa".0S33.  NODES"  1*  X.  X*X.  T • Y . X*X*Y*Y 


TABLE  13 


ARROftHEAO  *»()••  OELTA«.OSA.  MOOES*  !•  X*  X*X,  Y*t 


Table  15 

ftftftOVHEAO  *M»r  DfLT*".#**,  WOOES*  1 , Xt~4*X4  VYt 


AlUtOVHCAO  *•*>»«  K«l»|2t  K*|.t  DtLT*«.0*A.  MODES*  l»  X»  X*X 


0000*0  0000*0  0000*0  0000*0  0000*0 


•saooH  ,ctccc»o**o».\ai3o  • s**Un  •**•»  ovihoo 


I77|»OI  •2^047«0l  S*7)3SaOf  •7»4*«o»0*  1*8422*01  |*3»*7«01 


Table  21 


ARROWHEAD  A*R»  • M*I.S*2|,  K*l  • * DELTA* *05207 , MODES*  It 


1 •*ORO*OR  *S»S7S7*0R  -A.A2I3-0*  1«S*I*"0I  »#*MS«02 


Fig.  13  - Mesh  Used  for  Arrowhead  A = 4 
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